{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import cv2\n",
    "from scipy.interpolate import splprep, splev\n",
    "%matplotlib inline\n",
    "\n",
    "import triangle\n",
    "import triangle.plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## load image, find boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the image and convert to B&W\n",
    "I = plt.imread('hawk.png')\n",
    "I = I.sum(axis=2)\n",
    "I = np.where(I!=0, 1, 0)\n",
    "I = np.uint8(I)\n",
    "\n",
    "xmax = I.shape[1]\n",
    "ymax = I.shape[0]\n",
    "\n",
    "# blur, exand, then contract the image to smooth the lines\n",
    "blur=((3,3),1)\n",
    "erode=(3, 3)\n",
    "dilate=(5, 5)\n",
    "\n",
    "I = cv2.GaussianBlur(I , blur[0], blur[1])\n",
    "I = cv2.erode(I, np.ones(erode))\n",
    "I = cv2.dilate(I, np.ones(dilate))\n",
    "I *= 255\n",
    "\n",
    "# find the contours\n",
    "I, contours, _ = cv2.findContours(I, cv2.RETR_LIST, cv2.CHAIN_APPROX_NONE)\n",
    "\n",
    "I2 = 0 * I.copy()\n",
    "contour = contours[0]\n",
    "for contour in contours:\n",
    "    for i, j in contour[:,0]:\n",
    "        I2[j,i] = 127\n",
    "\n",
    "# get the image as coordinates\n",
    "tmpy, tmpx = np.where(I>0)\n",
    "IXY = np.vstack((tmpx, tmpy)).T\n",
    "IXY = np.array([0,1]) - IXY / np.array([-xmax, ymax])\n",
    "\n",
    "#only 4 matters (the longest)\n",
    "#rescale XY\n",
    "XY = np.array([0,1]) - contours[4][:,0] / np.array([-xmax, ymax])\n",
    "#XY = XY[:-1,:]\n",
    "\n",
    "# make a spline\n",
    "n = 500\n",
    "tck, u = splprep(XY.T, u=None, s=0.0)\n",
    "u_new = np.linspace(u.min(), u.max(), n)\n",
    "x_new, y_new = splev(u_new, tck, der=0)\n",
    "\n",
    "XYnew = np.vstack((x_new, y_new)).T\n",
    "nlist = np.arange(0,XYnew.shape[0])\n",
    "segments = np.vstack((nlist, np.roll(nlist, -1))).T\n",
    "\n",
    "#plt.scatter(IXY[:,0], IXY[:,1])\n",
    "plt.plot(XY[:,0], XY[:,1], 'r-', lw=4)\n",
    "plt.plot(x_new, y_new, 'b--s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make triangulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = triangle.triangulate({'vertices': XYnew, 'segments': segments}, 'pq30a.01D')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "triangle.plot.plot(plt.axes(), **t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Try PyAMG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.sparse as sparse\n",
    "triangles = t['triangles']\n",
    "\n",
    "row = triangles.ravel()\n",
    "col = np.roll(triangles, -1, axis=1).ravel()\n",
    "data = 1.0*np.ones((len(row),))\n",
    "\n",
    "L = sparse.coo_matrix((data, (row, col)))\n",
    "# symmetrize\n",
    "L = L.T + L\n",
    "# set to all minus ones\n",
    "L.data[:] = -1\n",
    "# zero diagonal\n",
    "L.setdiag(0.0)\n",
    "# collect totals\n",
    "diag = -np.array(L.sum(axis=1).ravel()).ravel()\n",
    "# set to minus totals\n",
    "L.setdiag(diag)\n",
    "# find Dinv\n",
    "n = L.shape[0]\n",
    "Dinv = sparse.spdiags(1.0/diag, [0], n, n, format='csr')\n",
    "L = Dinv.dot(L.dot(Dinv))\n",
    "L = L.tocsr()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyamg\n",
    "ml = pyamg.smoothed_aggregation_solver(L, max_levels=2, keep=True, max_coarse=10)\n",
    "\n",
    "AggOp = ml.levels[0].AggOp.tocsc()\n",
    "edges = []\n",
    "for c in range(AggOp.shape[1]):\n",
    "    row = AggOp.getcol(c).indices\n",
    "    # add to edges any edge in the aggregate\n",
    "    for i in row:\n",
    "        for j in row:\n",
    "            if i != j and L[i,j]:\n",
    "                edges.append([i,j])\n",
    "\n",
    "#AggOp = AggOp.T                \n",
    "AggOp = AggOp.tocsr()\n",
    "col = AggOp.indices\n",
    "aggtris = [t for t in triangles if col[t[0]]==col[t[1]] and col[t[1]]==col[t[2]]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(triangles.shape)\n",
    "print(len(aggtris))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make logo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = t['vertices'][:,0]\n",
    "y = t['vertices'][:,1]\n",
    "triangles = t['triangles']\n",
    "\n",
    "zfaces = np.array([np.linalg.norm([x[t[0]], y[t[0]]]) for t in triangles])\n",
    "plt.figure(dpi=200)\n",
    "plt.tripcolor(x, y, triangles, facecolors=zfaces, edgecolors='b', cmap=plt.cm.binary)\n",
    "\n",
    "# plot aggregate edges\n",
    "#for e in edges:\n",
    "#    plt.plot([x[e[0]], x[e[1]]], [y[e[0]], y[e[1]]], 'r-', lw=1.0)\n",
    "   \n",
    "# plot aggregate triangles\n",
    "#plt.plot(x[triangles[1]], y[triangles[1]], 'bo')\n",
    "plt.tripcolor(x, y, aggtris, facecolors=zfaces.max()*np.ones((len(aggtris,))), cmap=plt.cm.Blues_r)\n",
    "\n",
    "# plot title\n",
    "plt.text(0.4, -0.1,\n",
    "         'RAPtor: parallel algebraic multigrid',\n",
    "         fontsize=28,\n",
    "         family='sans-serif',\n",
    "         style='italic',\n",
    "         color=[0.3, 0.3, 0.3],\n",
    "        fontweight='bold')\n",
    "plt.axis('equal')\n",
    "plt.axis('off')\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.savefig('raptor-logo.png',\n",
    "            dpi=200,\n",
    "            bbox_inches='tight',\n",
    "            pad_inches=0,\n",
    "            transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
